Fit tweedie/delta models to biomass density and extract AIC

Author

Max Lindmark

Published

September 9, 2023

Intro

Fit Tweedie or delta models models to biomass density of cod, flounder, plaice and dab (juveniles and adults) between 1993-2020 in the Baltic Sea using sdmMTB fit different oxygen and temperature-related covariates, compare AIC. Select the “best” covariate for further trend and velocity analysis.

Load packages & source functions

# Load libraries, install if needed
pkgs <- c("tidyverse", "readxl", "tidylog", "RCurl", "devtools",
          "kableExtra", "viridis", "RColorBrewer", "here", "sdmTMBextra") 

if(length(setdiff(pkgs,rownames(installed.packages()))) > 0){

    install.packages(setdiff(pkgs, rownames(installed.packages())), dependencies = T)
  
  }

invisible(lapply(pkgs, library, character.only = T))

# Packages not on CRAN or dev version
# remotes::install_github("pbs-assess/sdmTMB", dependencies = TRUE)
library(sdmTMB)

# Source code for map plots
# You need: # devtools::install_github("seananderson/ggsidekick") # not on CRAN; library(ggsidekick)
devtools::source_url("https://raw.githubusercontent.com/maxlindmark/pred-prey-overlap/main/R/functions/map-plot.R")
options(ggplot2.continuous.colour = "viridis")

# Set path
home <- here::here()

Read data

# Read data
d <- #readr::read_csv("https://raw.githubusercontent.com/maxlindmark/spatial-metabolic-index/main/data/clean/catch_clean.csv") |>
  readr::read_csv(paste0(home, "/data/clean/catch_clean.csv")) |>
  rename(X = x, Y = y) |> 
  pivot_longer(c(cod_adult, cod_juvenile, dab_adult, dab_juvenile, flounder_adult,
                 flounder_juvenile, plaice_adult, plaice_juvenile),
               names_to = "group", values_to = "density") |> 
  separate(group, into = c("species", "life_stage"), remove = FALSE) |> 
  drop_na(depth, temp, oxy, sal, density)
Rows: 12385 Columns: 26
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr  (5): country, haul_id, ices_rect, month_year, substrate
dbl (21): year, lat, lon, quarter, month, sub_div, oxy, temp, sal, depth, x,...

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
rename: renamed 2 variables (X, Y)

pivot_longer: reorganized (cod_adult, cod_juvenile, dab_adult, dab_juvenile, plaice_adult, …) into (group, density) [was 12385x26, now 99080x20]

drop_na: removed 9,570 rows (10%), 89,510 rows remaining
# Read metabolic parameter estimates and left_join
mi_pars <- #readr::read_csv("https://raw.githubusercontent.com/maxlindmark/spatial-metabolic-index/main/data/clean/mi_params.csv") |>
  readr::read_csv(paste0(home, "/data/clean/mi_params.csv")) |> 
  dplyr::select(n_o2, E_o2, A0_o2, species) #  TODO: remove the extra columns for pressure based parameters
Rows: 4 Columns: 10
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (1): species
dbl (9): temp, po2, o2, A0_o2, A0_po2, n_po2, E_po2, n_o2, E_o2

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
# TODO: for now we'll use plaice parameters for flounder, see "00_estimate_mi_params.Rmd"
mi_pars
# A tibble: 4 × 4
    n_o2  E_o2  A0_o2 species 
   <dbl> <dbl>  <dbl> <chr>   
1 -0.300 0.294 0.0451 plaice  
2 -0.300 0.294 0.0249 cod     
3 -0.300 0.294 0.0371 dab     
4 -0.300 0.294 0.0120 flounder
mi_pars <- mi_pars |> 
  mutate(A0_o2 = ifelse(species == "flounder",
                        filter(mi_pars, species == "plaice")$A0_o2,
                        A0_o2))
filter: removed 3 rows (75%), one row remaining
mutate: changed one value (25%) of 'A0_o2' (0 new NA)
d <- left_join(d, mi_pars, by = "species")
left_join: added 3 columns (n_o2, E_o2, A0_o2)
           > rows only in x        0
           > rows only in y  (     0)
           > matched rows     89,510
           >                 ========
           > rows total       89,510
# Read size csv to calculate the metabolic index
sizes <- #readr::read_csv("https://raw.githubusercontent.com/maxlindmark/spatial-metabolic-index/main/data/clean/sizes.csv") |> 
  readr::read_csv(paste0(home, "/data/clean/sizes.csv")) |>
  mutate(group = paste(species, name, sep = "_")) |> 
  dplyr::select(group, B)
Rows: 8 Columns: 5
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (2): species, name
dbl (3): mat_w, max_w, B

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
mutate: new variable 'group' (character) with 8 unique values and 0% NA
d <- left_join(d, sizes, by = "group")
left_join: added one column (B)
           > rows only in x        0
           > rows only in y  (     0)
           > matched rows     89,510
           >                 ========
           > rows total       89,510
# Drop dab!
# d |>
#   filter(group == "dab_juvenile" & quarter == 1) |>
#   mutate(pres = ifelse(density > 0, "1", "0")) |>
#   ggplot(aes(X, Y, color = pres)) +
#   geom_point(size = 0.3) +
#   coord_fixed() +
#   facet_wrap(~year)
# 
# d |>
#   filter(group == "dab_juvenile" & quarter == 4) |>
#   mutate(pres = ifelse(density > 0, "1", "0")) |>
#   ggplot(aes(X, Y, color = pres)) +
#   geom_point(size = 0.3) +
#   coord_fixed() +
#   facet_wrap(~year)

d <- d |> filter(!species == "dab")
filter: removed 23,124 rows (26%), 66,386 rows remaining

Calculate the metabolic index

# Oxygen is ml/L, We want micro mol/L. 1 ml/l = 10^3/22.391 = 44.661 micro mol/l
d$oxy_si <- (d$oxy * (10^3)) / 22.391

# Calculate metabolic indices for a given mass and the temperature and oxygen in data
# Line 123 in https://github.com/fate-spatialindicators/SDM_O2/blob/master/code/mi_functions.R
kb <- 0.000086173324 # Boltzmann's constant
t_ref <- 15 # arbitrary reference temperature

# Calculate the metabolic index
d <- d |>
  mutate(inv_temp = (1/(temp + 273.15) - 1/(t_ref + 273.15)),
         phi = A0_o2*(B^n_o2)*oxy_si * exp((E_o2/kb)*inv_temp)) |> 
  drop_na(phi)
mutate: new variable 'inv_temp' (double) with 12,194 unique values and 0% NA
        new variable 'phi' (double) with 65,602 unique values and 0% NA
drop_na: no rows removed

Scale variables

d <- d |> 
  group_by(group) |> 
  mutate(phi_sc = as.vector(scale(phi)),
         temp_sc = as.vector(scale(temp)),
         temp_sq = temp_sc^2,
         oxy_sc = as.vector(scale(oxy)),
         oxy_sq = oxy_sc^2,
         sal_sc = as.vector(scale(sal)),
         depth_sc = as.vector(scale(depth)),
         depth_sq = depth_sc*depth_sc) |> # not sure this is needed!
  mutate(quarter_f = as.factor(quarter),
         year_f = as.factor(year)) |> 
  ungroup()
group_by: one grouping variable (group)
mutate (grouped): new variable 'phi_sc' (double) with 65,602 unique values and 0% NA
                  new variable 'temp_sc' (double) with 65,576 unique values and 0% NA
                  new variable 'temp_sq' (double) with 65,576 unique values and 0% NA
                  new variable 'oxy_sc' (double) with 65,568 unique values and 0% NA
                  new variable 'oxy_sq' (double) with 65,568 unique values and 0% NA
                  new variable 'sal_sc' (double) with 65,586 unique values and 0% NA
                  new variable 'depth_sc' (double) with 25,364 unique values and 0% NA
                  new variable 'depth_sq' (double) with 25,364 unique values and 0% NA
mutate (grouped): new variable 'quarter_f' (factor) with 2 unique values and 0% NA
                  new variable 'year_f' (factor) with 29 unique values and 0% NA
ungroup: no grouping variables
# Seems like a non-linear relationship with temperature, but a quadratic largely does the job??
ggplot(d, aes(temp_sc, density)) +
  geom_smooth(method = "lm", color = "grey50", se = FALSE) +
  geom_smooth(method = "gam", formula = y~s(x, k=4), color = "steelblue3", se = FALSE) +
  geom_smooth(method = "lm", formula = y ~ x + I(x^2), color = "tomato3", se = FALSE) +
  facet_wrap(~group, scales = "free")
`geom_smooth()` using formula = 'y ~ x'

ggplot(d, aes(oxy_sc, density)) +
  geom_smooth(method = "lm", color = "grey50", se = FALSE) +
  geom_smooth(method = "gam", formula = y~s(x, k=4), color = "steelblue3", se = FALSE) +
  geom_smooth(method = "lm", formula = y ~ x + I(x^2), color = "tomato3", se = FALSE) +
  facet_wrap(~group, scales = "free")
`geom_smooth()` using formula = 'y ~ x'

ggplot(d, aes(phi_sc, log(density + 1))) +
  geom_smooth(method = "lm", color = "grey50", se = FALSE) +
  geom_smooth(method = "gam", formula = y~s(x, k=4), color = "steelblue3", se = FALSE) +
  geom_smooth(method = "lm", formula = y ~ x + I(x^2), color = "tomato3", se = FALSE) +
  facet_wrap(~group, scales = "free")
`geom_smooth()` using formula = 'y ~ x'

# Plot mi by species and life stage
pal <- brewer.pal(name = "Dark2", n = 8)[c(7, 5, 3)]

d |> 
  mutate(species = str_to_title(species),
         life_stage = str_to_title(life_stage)) |> 
  ggplot(aes(phi, fill = species, color = species)) + 
  geom_density(alpha = 0.2) + 
  geom_vline(xintercept = 0, linetype = 2) + 
  coord_cartesian(expand = 0) + 
  facet_wrap(~ life_stage, ncol = 1) +
  theme_sleek(base_size = 9) + 
  theme(legend.position = c(0.85, 0.85)) +
  labs(x = "Metabolic index (\u03C6)", y = "Density", fill = "Species", color = "Species") + 
  scale_color_manual(values = rev(pal), name = "") +
  scale_fill_manual(values = rev(pal), name = "")
mutate: changed 66,386 values (100%) of 'species' (0 new NA)
        changed 66,386 values (100%) of 'life_stage' (0 new NA)

ggsave(paste0(home, "/figures/supp/mi_histogram.pdf"), width = 11, height = 11, units = "cm", device = cairo_pdf)

Fit models and save AIC

Fit models!

data_list_aic <- list()

for(i in unique(d$group)) {  
    
    dd <- d |> filter(group == i)
    
    mesh <- make_mesh(dd, xy_cols = c("X", "Y"), cutoff = 15) # Need specific meshes because not for all hauls did we record all species
    
    print(
    ggplot() +
      inlabru::gg(mesh$mesh) +
      coord_fixed() +
      geom_point(aes(X, Y), data = dd, alpha = 0.2, size = 0.5) +
      annotate("text", -Inf, Inf, label = paste("n knots = ", mesh$mesh$n), hjust = -0.1, vjust = 2) + 
      labs(x = "Easting (km)", y = "Northing (km)", title = str_to_sentence(str_replace(i, "_", " ")))
    )
    
    print(i)
    
    # 0. Linear oxygen and temperature, spatially varying quarter
    m0 <- sdmTMB(density ~ 0 + year_f + quarter_f + sal_sc + depth_sc + depth_sq + temp_sc + oxy_sc,
                 data = dd,
                 mesh = mesh,
                 family = tweedie(link = "log"),
                 spatiotemporal = "IID",
                 # this setting for spatial field(s) works better for majority of groups compare to single sf
                 spatial = "off",
                 spatial_varying = ~0 + quarter_f, # if enabled, keep quarter_f in main formula; 
                 time = "year")
    print("m0")
    sanity(m0)
    dd$m0_res <- residuals(m0)
    
    # 1. Linear oxygen and squared temperature, spatially varying quarter
    m1 <- sdmTMB(density ~ 0 + year_f + quarter_f + sal_sc + depth_sc + depth_sq + temp_sc + temp_sq + oxy_sc,
                 data = dd,
                 mesh = mesh,
                 family = tweedie(link = "log"),
                 spatiotemporal = "IID",
                 spatial = "off",
                 spatial_varying = ~0 + quarter_f,
                 time = "year")
    print("m1")
    sanity(m1)
    dd$m1_res <- residuals(m1)
    
    # 2. interaction between linear oxygen and temperature, drop quadratic
    m2 <- sdmTMB(density ~ 0 + year_f + quarter_f + sal_sc + depth_sc + depth_sq + temp_sc * oxy_sc,
                  data = dd,
                  mesh = mesh,
                  family = tweedie(link = "log"),
                  spatiotemporal = "IID",
                  spatial = "off",
                  spatial_varying = ~0 + quarter_f,
                  time = "year")
    print("m2")
    sanity(m2)
    dd$m2_res <- residuals(m2)

    # 3. breakpoint oxygen
    m3 <- sdmTMB(density ~ 0 + year_f + quarter_f + sal_sc + depth_sc + depth_sq + temp_sc + temp_sq + breakpt(oxy_sc),
                 data = dd,
                 mesh = mesh,
                 family = tweedie(link = "log"),
                 spatiotemporal = "IID",
                 spatial = "off",
                 spatial_varying = ~0 + quarter_f,
                 time = "year")
    print("m3")
    sanity(m3)
    dd$m3_res <- residuals(m3)

    # 4. logistic oxygen
    # m4 <- sdmTMB(density ~ 0 + year_f + quarter_f + sal_sc + depth_sc + depth_sq + temp_sc + temp_sq + logistic(oxy_sc),
    #              data = dd,
    #              mesh = mesh,
    #              family = delta_gamma(link1 = "logit", link2 = "log"),
    #              spatiotemporal = "off",
    #              spatial = "on",
    #              time = "year",
    #              control = sdmTMBcontrol(newton_loops = 2))
    # print("m4")
    # sanity(m4)
    # dd$m4_res <- residuals(m4)

    # 5. linear metabolic index
    m5 <- sdmTMB(density ~ 0 + year_f + quarter_f + sal_sc + depth_sc + depth_sq + phi_sc,
                 data = dd,
                 mesh = mesh,
                 family = tweedie(link = "log"),
                 spatiotemporal = "IID",
                 spatial = "off",
                 spatial_varying = ~0 + quarter_f,
                 time = "year")
    print("m5")
    sanity(m5)
    dd$m5_res <- residuals(m5)

    # 6. breakpoint metabolic index!
    m6 <- sdmTMB(density ~ 0 + year_f + quarter_f + sal_sc + depth_sc + depth_sq + breakpt(phi_sc),
                 data = dd,
                 mesh = mesh,
                 family = tweedie(link = "log"),
                 spatiotemporal = "IID",
                 spatial = "off",
                 spatial_varying = ~0 + quarter_f,
                 time = "year")
    print("m6")
    sanity(m6)
    dd$m6_res <- residuals(m6)
    
    # logistic metabolic index!
    # m7 <- update(m1, density ~ 0 + year_f + quarter_f + sal_sc + depth_sc + depth_sq + logistic(phi_sc))
    # print("m7")
    # sanity(m7)
    # dd$m7_res <- residuals(m7)
    
    # Plot residuals
    p1 <- dd |>
      dplyr::select(m0_res,
                    m1_res,
                    m2_res,
                    m3_res,
                    #m4_res,
                    m5_res,
                    m6_res,
                    #m7_res
                    ) |>
      pivot_longer(everything()) |>
      ggplot(aes(sample = value)) +
      stat_qq(size = 0.75, shape = 21, fill = NA) +
      facet_wrap(~name) +
      stat_qq_line() +
      labs(y = "Sample Quantiles", x = "Theoretical Quantiles") +
      theme(aspect.ratio = 1)
    
    print(p1)

    # Save AIC
    data_list_aic[[i]] <- AIC(m0, 
                              m1,
                              m2,
                              m3,
                              #m4,
                              m5,
                              m6#,
                              #m7
                              ) |>
      tibble::rownames_to_column("model") |>
      mutate(group = i)
    
}
filter: removed 54,050 rows (81%), 12,336 rows remaining
[1] "cod_adult"
[1] "m0"
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large
[1] "m1"
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large
[1] "m2"
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large
[1] "m3"
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large
[1] "m5"
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large
[1] "m6"
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large
pivot_longer: reorganized (m0_res, m1_res, m2_res, m3_res, m5_res, …) into (name, value) [was 12336x6, now 74016x2]

mutate: new variable 'group' (character) with one unique value and 0% NA
filter: removed 54,050 rows (81%), 12,336 rows remaining

[1] "cod_juvenile"
[1] "m0"
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large
[1] "m1"
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large
[1] "m2"
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large
[1] "m3"
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large
[1] "m5"
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large
[1] "m6"
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large
pivot_longer: reorganized (m0_res, m1_res, m2_res, m3_res, m5_res, …) into (name, value) [was 12336x6, now 74016x2]

mutate: new variable 'group' (character) with one unique value and 0% NA
filter: removed 54,802 rows (83%), 11,584 rows remaining

[1] "plaice_adult"
[1] "m0"
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large
[1] "m1"
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large
[1] "m2"
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large
[1] "m3"
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large
[1] "m5"
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large
[1] "m6"
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large
pivot_longer: reorganized (m0_res, m1_res, m2_res, m3_res, m5_res, …) into (name, value) [was 11584x6, now 69504x2]

mutate: new variable 'group' (character) with one unique value and 0% NA
filter: removed 54,802 rows (83%), 11,584 rows remaining

[1] "plaice_juvenile"
[1] "m0"
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large
[1] "m1"
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large
[1] "m2"
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large
[1] "m3"
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large
[1] "m5"
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large
[1] "m6"
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large
pivot_longer: reorganized (m0_res, m1_res, m2_res, m3_res, m5_res, …) into (name, value) [was 11584x6, now 69504x2]

mutate: new variable 'group' (character) with one unique value and 0% NA
filter: removed 57,114 rows (86%), 9,272 rows remaining

[1] "flounder_adult"
[1] "m0"
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large
[1] "m1"
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large
[1] "m2"
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large
[1] "m3"
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large
[1] "m5"
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large
[1] "m6"
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large
pivot_longer: reorganized (m0_res, m1_res, m2_res, m3_res, m5_res, …) into (name, value) [was 9272x6, now 55632x2]

mutate: new variable 'group' (character) with one unique value and 0% NA
filter: removed 57,112 rows (86%), 9,274 rows remaining

[1] "flounder_juvenile"
[1] "m0"
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large
[1] "m1"
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large
[1] "m2"
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large
[1] "m3"
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large
[1] "m5"
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large
[1] "m6"
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large
pivot_longer: reorganized (m0_res, m1_res, m2_res, m3_res, m5_res, …) into (name, value) [was 9274x6, now 55644x2]

mutate: new variable 'group' (character) with one unique value and 0% NA

# Save aic as data frames
data_aic <- bind_rows(data_list_aic)

write_csv(data_aic, paste0(home, "/output/data_aic_01.csv"))

Plot meshes in a single plot using patchwork

g1 <- d |> filter(group == unique(d$group)[1])
filter: removed 54,050 rows (81%), 12,336 rows remaining
mesh_g1 <- make_mesh(g1, xy_cols = c("X", "Y"), cutoff = 15)

g2 <- d |> filter(group == unique(d$group)[5])
filter: removed 57,114 rows (86%), 9,272 rows remaining
mesh_g2 <- make_mesh(g2, xy_cols = c("X", "Y"), cutoff = 15)

g3 <- d |> filter(group == unique(d$group)[3])
filter: removed 54,802 rows (83%), 11,584 rows remaining
mesh_g3 <- make_mesh(g3, xy_cols = c("X", "Y"), cutoff = 15)

g4 <- d |> filter(group == unique(d$group)[2])
filter: removed 54,050 rows (81%), 12,336 rows remaining
mesh_g4 <- make_mesh(g4, xy_cols = c("X", "Y"), cutoff = 15)

g5 <- d |> filter(group == unique(d$group)[6])
filter: removed 57,112 rows (86%), 9,274 rows remaining
mesh_g5 <- make_mesh(g5, xy_cols = c("X", "Y"), cutoff = 15)

g6 <- d |> filter(group == unique(d$group)[4])
filter: removed 54,802 rows (83%), 11,584 rows remaining
mesh_g6 <- make_mesh(g6, xy_cols = c("X", "Y"), cutoff = 15)


p1 <- ggplot() +
  inlabru::gg(mesh_g1$mesh) +
  coord_fixed() +
  geom_point(aes(X, Y), data = g1, alpha = 0.2, size = 0.5) +
  annotate("text", -Inf, Inf, label = paste("n knots = ", mesh_g1$mesh$n), hjust = -0.1, vjust = 2, size = 2.5) + 
  labs(x = "Easting (km)", y = "Northing (km)", title = str_to_sentence(str_replace(g1$group[1], "_", " "))) +
  theme_sleek(base_size = 8) +
  theme(axis.title.x = element_blank(),
        plot.margin = unit(c(0,0,0,0), "cm"))
  
p2 <- ggplot() +
  inlabru::gg(mesh_g2$mesh) +
  coord_fixed() +
  geom_point(aes(X, Y), data = g2, alpha = 0.2, size = 0.5) +
  annotate("text", -Inf, Inf, label = paste("n knots = ", mesh_g2$mesh$n), hjust = -0.1, vjust = 2, size = 2.5) + 
  labs(x = "Easting (km)", y = "Northing (km)", title = str_to_sentence(str_replace(g2$group[], "_", " "))) +
  theme_sleek(base_size = 8) +
  theme(axis.title.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.title.y = element_blank(),
        plot.margin = unit(c(0,0,0,0), "cm"))

p3 <- ggplot() +
  inlabru::gg(mesh_g3$mesh) +
  coord_fixed() +
  geom_point(aes(X, Y), data = g3, alpha = 0.2, size = 0.5) +
  annotate("text", -Inf, Inf, label = paste("n knots = ", mesh_g3$mesh$n), hjust = -0.1, vjust = 2, size = 2.5) + 
  labs(x = "Easting (km)", y = "Northing (km)", title = str_to_sentence(str_replace(g3$group[], "_", " "))) +
  theme_sleek(base_size = 8) +
  theme(axis.title.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.title.y = element_blank(),
        plot.margin = unit(c(0,0,0,0), "cm"))

p4 <- ggplot() +
  inlabru::gg(mesh_g4$mesh) +
  coord_fixed() +
  geom_point(aes(X, Y), data = g4, alpha = 0.2, size = 0.5) +
  annotate("text", -Inf, Inf, label = paste("n knots = ", mesh_g4$mesh$n), hjust = -0.1, vjust = 2, size = 2.5) + 
  labs(x = "Easting (km)", y = "Northing (km)", title = str_to_sentence(str_replace(g4$group[], "_", " "))) +
  theme_sleek(base_size = 8) +
  theme(axis.title.x = element_blank(),
        plot.margin = unit(c(0,0,0,0), "cm"))

p5 <- ggplot() +
  inlabru::gg(mesh_g5$mesh) +
  coord_fixed() +
  geom_point(aes(X, Y), data = g5, alpha = 0.2, size = 0.5) +
  annotate("text", -Inf, Inf, label = paste("n knots = ", mesh_g5$mesh$n), hjust = -0.1, vjust = 2, size = 2.5) + 
  labs(x = "Easting (km)", y = "Northing (km)", title = str_to_sentence(str_replace(g5$group[], "_", " "))) +
  theme_sleek(base_size = 8) +
  theme(axis.title.y = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        plot.margin = unit(c(0,0,0,0), "cm"))

p6 <- ggplot() +
  inlabru::gg(mesh_g6$mesh) +
  coord_fixed() +
  geom_point(aes(X, Y), data = g6, alpha = 0.2, size = 0.5) +
  annotate("text", -Inf, Inf, label = paste("n knots = ", mesh_g6$mesh$n), hjust = -0.1, vjust = 2, size = 2.5) + 
  labs(x = "Easting (km)", y = "Northing (km)", title = str_to_sentence(str_replace(g6$group[], "_", " "))) +
  theme_sleek(base_size = 8) +
  theme(axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        plot.margin = unit(c(0,0,0,0), "cm"))


p1 + p2 + p3 + p4 + p5 + p6

ggsave(paste0(home, "/figures/supp/mesh.pdf"), width = 17, height = 11, units = "cm", device = cairo_pdf)